library(tidyverse)## ── Attaching packages ──────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.0 ✓ purrr 0.3.3
## ✓ tibble 2.1.3 ✓ dplyr 0.8.4
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0## ── Conflicts ─────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(ggplot2)
library(plotly)## Warning: package 'plotly' was built under R version 3.6.2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutGet bird fragmentation data
bird_frag <- read.csv(file = "https://mlammens.github.io/ENS-623-Research-Stats/data/Logan_Examples/Chapter9/Data/loyn.csv")
summary(bird_frag)## ABUND AREA YR.ISOL DIST
## Min. : 1.50 Min. : 0.10 Min. :1890 Min. : 26.0
## 1st Qu.:12.40 1st Qu.: 2.00 1st Qu.:1928 1st Qu.: 93.0
## Median :21.05 Median : 7.50 Median :1962 Median : 234.0
## Mean :19.51 Mean : 69.27 Mean :1950 Mean : 240.4
## 3rd Qu.:28.30 3rd Qu.: 29.75 3rd Qu.:1966 3rd Qu.: 333.2
## Max. :39.60 Max. :1771.00 Max. :1976 Max. :1427.0
## LDIST GRAZE ALT
## Min. : 26.0 Min. :1.000 Min. : 60.0
## 1st Qu.: 158.2 1st Qu.:2.000 1st Qu.:120.0
## Median : 338.5 Median :3.000 Median :140.0
## Mean : 733.3 Mean :2.982 Mean :146.2
## 3rd Qu.: 913.8 3rd Qu.:4.000 3rd Qu.:182.5
## Max. :4426.0 Max. :5.000 Max. :260.0Log transform AREA
bird_frag$LOG_AREA <- log(bird_frag$AREA)x_values <- bird_frag$LOG_AREA %>%
round(2)
y_values <- bird_frag$GRAZE %>%
round(2)
z_values <- bird_frag$ABUND %>%
round(2)# Define regression plane -------------------------------------------------
# Construct x and y grid elements
x_grid <- seq(from = min(bird_frag$LOG_AREA), to = max(bird_frag$LOG_AREA), length = 50)
y_grid <- seq(from = min(bird_frag$GRAZE), to = max(bird_frag$GRAZE), length = 50)# Construct z grid by computing
# 1) fitted beta coefficients
# 2) fitted values of outer product of x_grid and y_grid
# 3) extracting z_grid (matrix needs to be of specific dimensions)
beta_hat <- bird_frag %>%
#lm(ABUND ~ LOG_AREA + GRAZE, data = .) %>%
lm(ABUND ~ LOG_AREA + GRAZE + log10(DIST) + log10(LDIST) + YR.ISOL + ALT, data = .) %>%
coef()
fitted_values <- crossing(y_grid, x_grid) %>%
#mutate(z_grid = beta_hat[1] + beta_hat[2]*x_grid + beta_hat[3]*y_grid)
mutate(z_grid = beta_hat[1] + beta_hat[2]*x_grid + beta_hat[3]*y_grid +
beta_hat[4]*mean(log10(bird_frag$DIST)) + beta_hat[5]*mean(log10(bird_frag$LDIST)) +
beta_hat[6]*mean(bird_frag$YR.ISOL) + beta_hat[7]*mean(bird_frag$ALT))
z_grid <- fitted_values %>%
pull(z_grid) %>%
matrix(nrow = length(x_grid)) %>%
t()# Define text element for each point in plane
text_grid <- fitted_values %>%
pull(z_grid) %>%
round(3) %>%
as.character() %>%
paste("ABUND: ", ., sep = "") %>%
matrix(nrow = length(x_grid)) %>%
t()# Plot using plotly -------------------------------------------------------
plot_ly() %>%
# 3D scatterplot:
add_markers(
x = x_values,
y = y_values,
z = z_values,
marker = list(size = 5),
hoverinfo = 'text',
text = ~paste(
"ABUND: ", z_values, "<br>",
"GRAZE: ", y_values, "<br>",
"LOG_AREA: ", x_values
)
) %>%
# Regression plane:
add_surface(
x = x_grid,
y = y_grid,
z = z_grid,
hoverinfo = 'text',
text = text_grid
) %>%
# Axes labels and title:
layout(
title = "3D scatterplot and regression plane",
scene = list(
zaxis = list(title = "y: ABUND"),
yaxis = list(title = "x2: GRAZE"),
xaxis = list(title = "x1: LOG_AREA")
)
)